'''
Purpose: To detrend a matrix of time series (inputMatrix) using 
regression of against a polynomial time trend of order p

Usage: residual = detrend(inputMatrix, p)
if p == -1    return the input
if p == 0    subtract the mean
if p > 0     detrend

Created on Aug 16, 2010

@author: Peter Harrington
peter.b.harrington@gmail.com
'''
from numpy import *
def detrend(inputMatrix, p):
    if p == -1:
        return inputMatrix
    
    numObs = inputMatrix.shape[0]
    u = matrix(ones((numObs,1)))                #create 0th order term
    
    if (p > 0): #create a time matrix
        timep = matrix(zeros((numObs,p)))
        t = matrix(range(numObs))
        print t
        tp = t.transpose()/float(numObs)
        print tp
        m = 1
        while (m <= p):
            timep[:,m-1] = multiply(tp,m)
            m += 1
        print timep
        xmat = concatenate((u,timep),1)
    else:
        xmat = u
    print xmat
    xpxi = (xmat.transpose()*xmat).I            #denom for ols
    beta = xpxi*(xmat.transpose()*inputMatrix)  #calc ols
    residuals = inputMatrix - xmat*beta         #subtract residuals
    return residuals